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Abstract 

We present an approach for the description of fluctuations that are due to finite system size 
induced correlations in the Kuramoto model of coupled oscillators. We construct a hierarchy for 
the moments of the density of oscillators that is analogous to the BBGKY hierarchy in the kinetic 
theory of plasmas and gases. To calculate the lowest order system size effect, we truncate this 
hierarchy at second order and solve the resulting closed equations for the two-oscillator correlation 
function around the incoherent state. We use this correlation function to compute the fluctuations 
of the order parameter, including the effect of transients, and compare this computation with 
numerical simulations. 



1 



Systems of coup 



phenomena 



ed oscillators appear as models for the dynamics of a wide range of 



2, 3, 4], 5], Q]. The Kuramoto model is a simple and oft-studied description 



of coupled oscillators which, in the limit of an infinite number of oscillators, exhibits a 



phase transition from an incoherent state to phase locked dynamics , lUj . However, 

numerical simulations show the appearance of fluctuations that are due to finite system size 
effects even in the absence of any external noise. Because the system is deterministic, these 
fluctuations are a manifestation of multi-oscillator correlations and are expected to vanish 
in the infinite oscillator limit, with potentially divergent behavior near the transition jl^ . 
While there has been some effort towards an analytic treatment of the fluctuations in the 
Kuramoto model {13I . there is at present no systematic approach. Here, we present a 
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161 ]. Our methods 



statistical formalism which draws upon the kinetic theory of plasmas 
are generalizable to any oscillator model. 

The Kuramoto model describes the phase evolution of N oscillators and is given by 

K N 

^Wi + TuEm-^ i = l,...,N, (1) 

where K is the coupling strength; the u>i are drawn from a distribution g(u), assumed to 
be symmetric and of zero mean. The coupling function f(8) can be any function. In the 
original Kuramoto model f(8) = sin 9, which we use for our simulations. 



In the N — > 00 limit, Kuramoto showed |8j that as the coupling K is increased from 0, 
this model exhibits a phase transition described by the order parameter 

1 N 

Z = -Y e»i = re** (2) 

which is a measure of the level of synchrony in the population. Kuramoto found a continuous 
transition from a phase of complete incoherence (r = 0) in the population to a relative 
degree of coherence (r > 0) for K greater than K c = 2/7177(0). However, for a finite number 
of oscillators, r will fluctuate. Figure [1] shows (r 2 ) (where (•) represents an ensemble average 
over frequencies and initial angles) for a numerical simulation of N = 50 oscillators. We see 
that the fluctuations smooth the sharp transition from incoherence to coherence. 

As we will show, typical with phase transitions, the correlations become enhanced near 
the onset of the transition (critical point). At low K, (r 2 ) « 1/N, consistent with the finite 
size effects for the free (K = 0) model. It reasonable to suppose that in the incoherent state 
all correlation effects (when K 7^ K c ) are finite size effects due to the coupling, i.e. the 
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FIG. 1: Phase transition in the Kuramoto model from asynchronous (K < K c ) to synchronous 
[K > K c ) behavior. The solid line is the mean field prediction for r 2 ; crosses are simulation data. 

homogeneous all-to-all connectivity forces the suppression of fluctuations as N — > oo. One 
of our goals is to calculate (r 2 ) including the fluctuations due to finite N in the incoherent 
state. 

Mirollo and Strogatz [l^] analyzed the stability of the incoherent state using a Fokker- 
Planck formalism. In the absence of external additive noise, their Fokker-Planck equation 
has the form of a continuity equation. They found that the incoherent state has a continuum 
of marginally stable modes, which are made stable by additive noise. In the ensuing, we 
will generate a series of equations analogous to the BBGKY hierarchy for which the Mirollo- 
Strogatz continuity equation is the truncation at first order. Our strategy is to consider an 
expansion using 1/N as a small parameter. 

The complete oscillator probability density 

1 N 

n(9, u,t) = -'£ 5(9 - 9 l (t))5(uj - u t ) (3) 

^ i=l 

satisfies the continuity equation 

^ + = -K^ jT fj f{9> - 9)n(9', u/, t)n(9, u>, t)d9'dco'. (4) 

Equation (j4]) is analogous to the Klimontovich equation in the plasma context and is still 
an exact description of the microscopic dynamics. The complete oscillator distribution ([3]) 
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satisfies the Klimontovich equation (J4]) exactly if the Kuramoto system flTJ is obeyed. Solving 
the Klimontovich equation for the complete distribution is equivalent to solving the original 
system and is equally difficult. The strategy of kinetic theory is to consider the smoothed 
probability density functions of the oscillators by taking ensemble averages. 

The one-particle probability density function (PDF) (first moment of n(8,uj,t)) is given 

by 

p 1 (6,u,t) = (n(6,u,t)) (5) 

where brackets denote the ensemble average over initial conditions and frequencies. The 
density pidOdu represents the mean fraction of oscillators within frequency range (ou, ou + du) 
and angle range (6,6 + dd). We note that J 27r pi(6,uj,t)d6 = g(uj). Henceforth, we will use 
the compact notation x = (6,u). Taking the expectation value of Eq. (J4]) gives 

= -K— j ^ ^ f(6' - 6)C(x; x', t)d9 ' dJ (6) 

where 

C(x; x', t) = p 2 {x\ x', t) - pi(x, t)pi(x', t) (7) 
is the "connected" two oscillator correlation or moment function with 

(n(x, t)n(x', t)) = p 2 {x\ x', t) + j^d( x - x ')pi{x, ( 8 ) 

The self-fluctuation term drops out in Eq. because we consider f(0) = 0. 

The RHS of ([6]) describes two oscillator interactions and is comparable to the collision 
integral from the kinetic theory of gases and plasmas. Neglecting the collision integral leads 
to the Vlasov equation, which amounts to a mean field approximation. The Vlasov equation 
and corresponding Fokker-Planck equation, which includes a diffusive term when external 



noise is mc 
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uded, has been studied for coupled oscillators previously in many contexts 
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231 ] . Although the Vlasov equation has the same form as Eq. (j4j), 
the two should not be confused. p x (x, t) is a smooth function representing the expectation 
value of the number density over initial conditions and frequencies, whereas n(x, t) is an 
operator-valued distribution and contains all statistical information about the system. 

We obtain an equation for C(x;x',t) by multiplying Eq. (j4j) by n(x',t) and taking the 
expectation value. This will result in an equation that depends on the three oscillator mo- 



ment function. Continuing this process for higher moments results in the BBGKY hierar- 



chy 
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We truncate the hierarchy at second order, expecting the correlation C(x; x', t) 



to be Oil IN) and a general connected n-point function to be Oil/N n x ) as is consistent 

nn 

with previous simulations [12|, H4[ . 

Using Eq. §§§ and removing terms expected to be 0(1/N 2 ) yields 

4t + ^ mi + + K L Oik m ~ h)+ k m - <>^ c ^ *». «> 

r 2« Q 



/ " 9 l )p l (x 1 ,t)C(x 2 , x 3 , t)d9 3 du 3 

-oo JO OUi 



lo d0 1 

r 2n Q 



/OO Q 
/ Tuff^ - e 2 ) Pl (x 2 , t)C(x 3 , x 1: t)d6 3 doj 3 } 
-oo JO OVo 



= " fi) + " fc)M*i.*)Pi(*2,t), (9) 

Equations ([H]) and form a Gaussian closure of a kinetic theory describing the Kuramoto 
model. We use this to calculate the fluctuations about the incoherent state. We start with 
the ansatz 151 ]: 

/oo p2tt rt 

du)' Y duJ 2 \ dQ' x dQ' 2 \ dtx(x[, x' 2 , t') 
-oo J o Jt 

x P(x 1 ,x' 1 ,t-t')P(x 2 ,x 2 ,t-t'). (10) 

where the initial conditions are imposed at t and t < t' < t. Using Eq. ffTUj) in Eq. Q), we 
obtain the dynamics for the propagator P, 

d d „ d r°° r 2w 



arid f°° f zn 
{^7 + ^7^ + K T^ / / f(02-d 1 ) Pl (x 2 ,t)d9 2 dw 2 }P(x 1 ,x' 1 ,t-t') 
at oUi ov\ J~oo Jo 

r) roo p2tv 

+ K i^r / f(02-e 1 )p 1 (x lt t)p(x 2 ,x? 17 t-t)de 2 du 2 = o, (11) 

OV\ J-oo JO 



where 

X (x 1 ,x 2 ,t) = -^[A/(0 2 _ ft) + J_/(ft _ ft)]^,*)^,*) (12) 

and the initial condition is P(x; x', 0) = 5(x — x'). 

The fluctuations in the order parameter are given by 

(r 2 ) = (ZZ*) = f dudu'dOdO' Xcu,#,t)n(cy,#',t))e l(f? - 9,) (13) 



We consider fluctuations in the incoherent state and thus seek solutions to §§§ and 
such that pi(x, t) = ^-g{uj). From Eqs. ([7]) and (jSj), we see that a computation of the 
fluctuations amounts to a calculation of the connected correlation function, which is phase 



invariant (because p\ is independent of 9), so that C(9i,9 2 ,u>i,u> 2 ,t) = C{9\ — 9 2 ,ui,u 2 ,t). 
Hence, the collision integral in Eq. ([6]) is zero, making pi(9,uj) = g{uj)/2-K an exact solution 
of the equations. Taking the Fourier and Laplace transforms in 9 and time of Eq. (TTU|) gives 

tiK\th\ f 1 f°° f f 

C n (u u u 2 ,s) = — -— — / du[du' 2 dT ds! ds 2 g{u[)g(uj' 2 ) (14) 

Z7T iV J —oo J Ci J C2 

x P_ n (u h u[, Sl )P n (uj 2 , u>' 2 , s 2 )-e^ +s ^ s) \ 

s 

where n is the Fourier mode index, s 12 is a Laplace transform variable and r = t — t' . By 
definition, the Laplace contours C\ and C 2 are arranged such that they are to the right of 
all poles in P_ n and P n , respectively. Using Eqs. ([7]), ([8]) and ffT4"l) in Eq. f|T3|) gives 

(r 2 (r)) = 4tt 2 | dcudu'C^cu, u/, r) + (15) 

because (Z) = 0. 

We can obtain a general expression for (r 2 ) without explicitly solving for the correlation 
function. From equation ( TTTT) . we can derive the relation 

where 

A n s = 1 + in-K/ n / — — — 17 
J -oo s + 277,0; 

is the analog of the dielectric response function. Using Eqs. (fTEl) and (fl4l in Eq. (|T5|) yields 



{r [T)) iKNnJc Axis-so) K6S 



-e ST (18) 



Ai(s 

where so is the zero of A n (s). The strategy of the calculation leading to Eq. ( |T8l) is similar 



to the calculation of the Lenard-Balescu collision integral [15l. Il6l|. 

For the specific frequency distribution g{uj) = (7 / V) (l/(u; 2 + 7 2 )) (i.e. a Lorentz distri- 
bution), Eq. (|T8|) evaluates to 

{ [ )) N K c — K N K c — K 1 > 

where K c = 27 for the Lorentz frequency distribution. (r 2 (0)) = 1/N because the initial 
conditions for Eqs. ([6]) and (J9J) are such that pi (x, 0) is the equilibrium incoherent state and 
C(xi,x 2 ,0) = 0. For the uncoupled system K = 0, so (r 2 ) = 1/N as expected. We also see 
that the amplitude of the fluctuations and the transient decay time become singular at the 
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critical point K = K c . At criticality, we obtain the expression (r 2 (r)) = (l/iV)(l + K c r). 
The closer K is to criticality, the less this calculation should be valid. Near critical behavior 
requires an analysis of all orders in the 1/N expansion. Dynamically, the implication is that 
as the coupling strength nears criticality, oscillators will interact more strongly and higher 
order correlations will become more important. The result (r 2 (oo)) = K C /[N(K C — K)] was 



first derived by Daido [14J with a completely different approach. Our method facilitates 
a systematic expansion in 1/N, in addition to providing an examination of the transient 
behavior of (r 2 (r)). 

We can examine the transient behavior of the correlations by solving Eq. (Tl4l) for the 
Lorentz distribution. We first solve for the propagator in Eq. ( Hip by taking a Fourier series 
expansion in 9 and Laplace transform in time, to obtain 

PnM, S) = 1 6(Ul ~ U[) , ^P^^ „ (20) 

where s is the Laplace transform variable and A„(s) = 1 — (K/2)\n\/(s + \n\~f). The prop- 
agator (120]) has poles along the imaginary axis corresponding to the continuous spectrum 
of marginally stable modes as well as those given by the discrete zeros of A n (s), which for 



K < K c = 27 are real and negative 



We then use Eq. (|20|) in Eq. (Tl4|) and take the inverse Laplace transform. For the n = 1 
mode, this gives 

(iOJ + f)(-iuj' + fO / e (i(oj-a/)"r) x \ 



K 1 

Cl(^,^',r) = ——g(w)g(u)) 



+ 2(f-f + ^)((f -f) 2 + H 2 )^ e 

ZT 2 1 
+ _ e -(Kc-K)r_ 1 \ 



(21) 



4(X - 2f ) (f - f - iw)(f - f + io;/) 

where r = £ — t - The other modes are given by C_i = Cjf and C n = for n 7^ ±1 since 
/(#) = sin#. The correlation function will thus have the form 

C(u, J, 6 - 9', r) = djV-V + (22) 

Only the last term in Eq. (I2T1) contributes to the transient in Eq. (fT9"I) . The correlation 
function contains modes which consist of all possible pairings of a marginal oscillating mode 



with the decaying mode. While the marginal modes do not decay in the correlations, they 
do not affect the decay of (Z(t)) because of a Landau damping-like dephasing effect that is 
described in Ref. {24]. The marginal modes also have no effect upon the transient behav- 
ior of (r 2 (r)). We should expect a similar result for higher moments. At this order, the 
marginal modes are not rendered stable by finite size effects as they are with the addition 
of external additive noise 17j. Should stabilization occur due to the intrinsic fluctuations, 



it will necessarily be a consequence of higher order effects. 

We compare our analytical results to numerical simulations of the Kuramoto system. 
Figure [2] a) shows the asymptotic value of (r 2 ) for various values of K and N. The analytical 
prediction matches extremely well for N = 500 and reasonably well for N = 50 and N = 100. 
Only at N = 10 are there significant deviations from the prediction. Fig. [2] b) shows the 
transient behavior of (r 2 ). The results match quite well below K/K c = 0.8. Numerical 
results for the correlation function integrated over ui, ui' are shown in Fig. [3j The simulation 
agrees well with the prediction Eq. ( 122]) except near the critical point as expected. 

Our calculation is the first presentation of a systematic approach to understanding the 
fluctuations due to finite size effects to an arbitrary order in 1/N. Although we truncate 
at lowest order, our approach allows a truncation at any level of the moment hierarchy 
to produce an expansion in 1/N. We note that Ref. 25[ found that when the oscillators 
are driven with Gaussian noise, 1/N dependence is still seen in the fluctuations of the order 
parameter. Additionally, our methodology could be used to study the evolution of the phase 
of Z, as in Ref. 26]. 



Some previous work 
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25 



271 ] for both phase and pulse coupled oscillators 
also start with a continuity equation similar to Eq. (j4j) but either go directly to mean field 
theory, with and without an external noise source to approximate fluctuations, or assume 
the fluctuations are Gaussian. References 
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derive a kinetic theory for a network of 



integrate-and-fire neurons by constructing a moment hierarchy similar to ours that is closed 
using the maximum entropy principle. However, this work differs from ours in that the 
hierarchy is built from a Boltzmann-like equation for a one-particle distribution function 
with stochastic inputs and hence does not capture the same correlation effects that we find 
by starting from a continuity equation that contains the full statistics of the system. 

We feel it is important to stress that the Klimontovich continuity equation (Eq. 
is not an approximation. The approximation appears in the method of finding solutions. 
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FIG. 2: a) Simulated and predicted N(r 2 } vs. K/K c for various values of N for large times. The 
data are evolved for 6 time constants at each K (r = \/{K c — K)) and averaged over 1000 different 
initial conditions and frequencies. The frequency distribution is Lorentz with 7 = 0.05. The 
simulation was performed with a time step of 0.05. The initial distribution of angles was uniform, 
b) Time evolution of (r 2 (r)) vs. r for various values of K and for N = 100. At each time point the 
values are averaged over 10,000 different initial conditions and frequencies. All other parameters 
as above. 

The mean field limit is equivalent to setting correlations to zero. Computing the moment 
hierarchy allows for an expansion which accounts for the effects of correlations. We produce 
a systematic method for deriving such an expansion and show explicitly in what regime 
higher order correlations can be ignored. We also note that the kinetic theory approach 
is only one avenue to understanding correlations. An alternative formulation is through 
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FIG. 3: C(uj,lo', 9 — 9') integrated over oj and oj' versus 9 — 9' for N = 100 and various values of 
K. Frequency distribution is Lorentz with 7 = 0.05 and the initial angle distribution is uniform. 
Results are averaged over 100,000 samples in a 2 dimensional histogram with 100x100 bins. The 
data is then averaged over angle differences and then put into a one dimensional histogram with 
bins of width 10. The time step for the evolution is 0.05. 



a statistical field theory approach 28[], which facilitates the construction of an expansion 
without resorting to a moment hierarchy. 
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